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ABSTRACT  h yry  V/ORT) 


theoretical  analysis  wt'.s  developed  for  predicting  the  flow  properties 
in  the  mixing  region  between  a particle-laden,  turbulent  rocket  exhaust 
and  a surrounding  air  stream.  It  was  assumed  that  the  turbulent  boundary 
layer  equations,  modified  to  account  for  particles,  were  valid  within  the 
mixing  region.'-\In  treating  the  chemical  aspects  of  the  problem  it  was 
assumed  that  the  (flow  was  in  equilibrium  in  accordance  with  three  chemi- 
cal reaction  equations^'  The  chemical  species  comprising  the  flow  were 
limited  to  the  following:^  HgO,  Hg,  02,  CO,  C02,  N2,  HC1,  Al,  Al^,  and 
one  additional  inert  species  which  remained  as  an  £nput~s5lection^  The 
phenomenological  model  employed  for  the  turbulent  transport  coefficients 
is  discussed  in  detail  and  compared  with  various  other  models. 

The  solution  of  the  partial  differential  equations  was  obtained  by  trans- 
forming the  equations  using  the  von  Mises  transformation^) expressing  the 
transformed  equations  in  finite-difference  form,  and^eolving  the  resulting 
equations  utilizing  a computer  program  developed  for  the  SRU  1107. 


Results  from  the  computer  program  were  successfully  compared  with  experi- 
mental results  obtained  from  air-augmentation,  free  jet,  and  fuel  injection 
experiments.  Output  data  from  the  computer  program  comprises  velocity, 
temperature,  density,  species  concentration,  and  Mach  number  profiles 
at  various  axial  locations  along  the  mixing  region. 
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NOTATION 


g 

h 

A h* 

J 

K 

P 

kf 

L 

1 

Le 

• 

m 

“p 

M 

n 

P 

Pr 


9 

area,  ft 

width  of  mixing  zone,  ft. 
constant  in  eddy  viscosity  expression,  ft. 
concentration  of  species  i,  lbm  of  i/lbm  of  mixture 
specific  heat  at  constant  pressure,  B/ltom-'R 
molecular  diffusion  coefficient,  ft  /sec 
internal  energy,  B/lbm 

gravitational  constant,  32.17  lbm— ft/lbf— sec 
static  enthalpy,  B/lbm 
heat  of  formation,  B/lbm 
mechanical  equivalent  of  heat,  778  ft-lb/B 
equilibrium  constant 

forward  reaction  rate  constant,  ft^/lb-mole  sec 
reference  length,  ft. 

Pr&ndtl' 8 mixing  length  parameter,  ft. 
turbulent  Lewis  number 
mass  flux,  ltm/ft^-scc 
particle  density,  ltm  of  Pf ft^  of 
Mach  number 

exponent  in  eddy  viscosity  expression 
static  pressure,  psfa 
turbulent  Prandtl  number 

universal  gas  constant,  1544  ft-lb/lb-mole  *R 
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radius,  ft. 

turbulent  Schmidt  number 

static  temperature,  *R 

axial  velocity  component,  ft/sec 

radial  velocity  component,  ft/sec 

net  rate  of  production  of  species  i,  lbm/ft  -sec 

molecular  weight  of  species  i,  lbm/lb-mole 

axial  distance,  ft. 

radial  distance,  ft. 

exponent  for  two-dimensional  ( s =0)  or  axisynmetric 
flow  ( S =1) 

eddy  diffusivity,  ft^/sec 

gas  density,  lbm  of  gas/ftJ  of  gas 

particle  density,  lbm  of  p/t\r  of  gas 

similarity  parameter 

2 

shear  stress,  lb/ft 
stream  function,  lbm/sec 


Subscripts 


mass  diffusivity 


heat  diffusivity 


m,  n 


species 

grid  point  coordinates 
particle 

momentum  diffusivity 
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l.o  INTRODUCTION 


1.1  General  Remarks 

Wytlcal  stMiss  involving  the  turbulent  mixing  of  multi-component  gasee 
required  in  a number  of  flow  probleme.  Examplea  of  auch  problems  are 
those  pertaining  to  base  heating,  wakes  behind  bluff  bodies,  loss 
munication  du.  to  pl».  attenuation,  aM  most  recently,  air  augmentation. 
Fundamental  studies  of  the  miking  process  between  two  moving  streams  have 
been  cc Mooted  over  the  past  several  years  with  varying  degrees  of  success. 
However,  it  has  been  only  recently,  with  the  advent  of  computers, 

that  Significant  progress  has  been  achieved  in  analytically  p c ng 
properties  within  turbulent  mixing  regions.  An  excellent  mon  gr  p 
bribing  the  recent  theoretical  developments  pertaining  to  turbulent  let 
has  been  compiled  by  Abramovich  (D*. 

Recent  air-ausaentation  studies  indicate  that  accurate  prediction  of  engine 
performance  require,  a detailed  study  of  the  mixing  region  between  a rocket 

exhaust  jet  and  a confined  secondary  air  stream.  Since  rock 

•»N,oT  not  only  must  the  mixing  process 
gases,  in  general  contain  excess  fuel,  not 

be  analysed,  but  alec  chemical  reaction  muat  be  taken  into  consideration. 
The  problem  ia  further  complicated  by  the  fact  that  the  entire  flow  is 
confined  within  a duct,  thue  nec.seiteting  the  consideration  of  axial 
pressure  gradients.  Analyses  of  flows  similar  to  the.,  found  in  this 
type  of  problem  hsv.  rscsntly  h.sn  developed  by  Ubby  (a)  and  VaSiliu  O). 
Libby,  for  example,  ha.  analysed  the  unconfined  mixed  flew  field  for 
eupersonic  combustion  studies  easing  equilibria  chemistry*  whereas 


* ♦ „ w^fftrences  listed  in  the  Bibliography^ 

♦Numbers  in  parentheses  refer  to 
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Vesiliu  assumed  nonequilibrium  chemistry  in  his  analysis  o t the  mixing 
between  a rocket  exhaust  Jet  and  a supersonic  air  stream.  Mikhail  (4), 
on  the  other  hand,  has  treated  the  mixing  of  coaxial  incompressible  streams 
in  a duct  but  did  not  consider  chemical  reactions. 

1.2  Scope  of  Present  Analysis 

The  analysis  and  associated  computer  program  described  herein  were  developed 
specifically  for  the  air- augment  at  ion  problem.  However,  it  should  be  noted 
that  the  computer  program  may  be  employed  for  other  problems.  The  primary 
objective  of  the  analysis  was  to  provide  a method  for  calculating  the 
velocity,  temperature,  and  composition  profiles  in  the  mixing  region  between 
a fuel- rich  rocket  exhaust  and  an  air  stream  confined  within  a duct  of 
specified  geometry.  Since  the  exhaust  of  solid  propellant  rockets  often 
contains  solid  and/or  liquid  particles,  the  analysis  was  developed  for 
either  gas-particle  systems,  or  systems  consisting  solely  of  gaseous  species. 
In  treating  the  chemical  aspects  of  the  problem  it  is  assumed  that  the  flow 
is  in  chemical  equilibrium  iu  accordance  with  three  chemical  reaction 
equations;  or,  if  desired,  the  flow  may  be  treated  as  frozen.  The  total 
number  of  chemical  components  comprising  the  flow  is  limited  to  a maximum 
of  ten  species. 

The  flow  model  employed  for  the  analysis  is  presented  in  Figure  1.  At  the 
periphery  of  the  rocket  nozzle  exit,  where  the  air  stream  and  rocket 
exhaust  first  come  into  contact,  a mixing  region  is  assumed  to  begin.  As 
the  flow  moves  downstream,  the  mixing  region  widens  on  both  sides  of  the 
plume  slipline  - the  inner  boundary  approaching  the  duct  centerline  and  the 
outer  boundary  the  duct  wall.  At  3ome  distance  downstream  from  the  nozzle 
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the  mixing  region  may  be  considered  to  extend  across  the  entire  duct.  Iu 
establishing  the  coordinates  of  the  plume  slipline  the  method  of  character- 
istics solution  is  employed  (5).  In  the  areas  outside  the  mixing  region, 
the  flow  is  considered  to  be  inviscid  with  uniform  properties  in  the 
radial  direction.  Within  the  mixing  region  turbulent  transport  processes 
are  assumed  to  prevail  with  all  molecular  transfer  considered  negligible. 
This  assumption  is  justified  for  the  problem  under  consideration  since 
the  flow  is  generally  turbulent,  and  the  turbulent  transport  coefficients 
are  usually  at  least  an  order  of  magnitude  greater  than  their  molecular 
counterparts.  The  phenomenological  model  employed  for  the  turbulent  trans' 
port  coefficients  in  the  present  study  is  discussed  in  detail  and  compared 

with  various  other  models. 

The  major  assumptions  employed  in  the  development  of  the  analysis  are  as 
follows : 

1.  Hie  flow  is  either  axisymmetric  or  two-dimensional. 

2.  The  entire  flow  field  is  turbulent. 

3.  The  gas  obeys  the  perfect  gas  law. 

4.  The  radial  pressure  gradient  across  the  mixing  zone  is  negligible. 

5.  Mass  transfer  due  to  thermal  and  pressure  gradients  is  negligible. 

6.  The  turbulent  Prandtl  and  Lewis  numbers  of  the  gas  are  constant  but 


may  have  values  other  than  unity. 

7.  The  particles  may  have  a turbulent  Lewis  number  different  than  that 
of  the  gas. 

8.  The  velocity  and  temperature  lags  between  the  gas  and  particles 
are  negligible. 
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9.  The  volume  occupied  by  the  particles  is  negligible  compared  to  that 


of  the  gas. 

10.  The  boundary  layer  between  the  mixed  flow  and  the  duct  wall  is 
neglected.  Hie  boundary  layers  in  the  rocket  exhaust  and  tuc  air 
stream  at  the  initial  point  of  contact,  however,  may  be  considered. 

11.  The  turbulent  boundary  layer  equations  are  valid  within  the 
mixing  region. 

1.3  Computer  Program 

The  solution  of  the  partial  differential  equations  is  obtained  by  first 
transforming  the  equations  using  the  von  Mises  transformation,  expressing 
the  transformed  equations  in  finite-difference  form,  and  then  solving  the 
resulting  finite-difference  equations  utilizing  a computer  program 
developed  for  the  SKU  1107*  Details  of  the  computer  program  are  presented 
in  Reference  6.  Consideration  is  given  to  the  stability  and  convergence 
of  the  solution  utilizing  the  criteria  established  by  Wu  (7). 

The  parameters  which  must  be  input  to  the  program  are  as  follows:  (a)  the 

initial  conditions  of  velocity,  temperature,  pressure,  and  species  con- 
centrations of  the  two  streams,  (b)  axial  step-size,  (c)  duct  geometry, 

(d)  plume  slip-line  coordinates,  (e)  the  turbulent  Lewis  and  Prandtl 
numbers,  and  (f)  computer  control  information  as  noted  in  Reference  6. 

The  output  data  comprises  velocity,  temperature,  density,  species  concentra- 
tion, and  Mach  number  profiles  within  the  mixing  region  at  various  axial 
locations;  axial  pressure  distributions;  and  axial  thrust  acting  on  the 
duct. 
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The  major  limitations  of  the  program  are  &3  follows:  (a)  the  calculations 

are  valid  o"ly  in  the  temperature  range  ^50R  to  995C®>  (b)  the  gaseous 
species  which  may  considered  are  limited  to  CO,  C02>  &>»  ^ 2> 

Al;  in  addition,  provision  for  one  additional  inert  specie  has  been 
incorporated  in  the  program,  and  (c)  the  particle  specie  which  may  be 
considered  is  Al^O^. 

Typical  run  time  of  the  program  is  10  minutes  for  a case  involving  a rocket 
nozzle  radius  of  0,21  ft,  a duct  radius  of  0,56  ft,  axial  step-3ize  of 
0.01  ft,  and  a total  axial  distance  of  6.0  ft. 
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DERIVATION  OF  THE  GOVERNING  EQUATIONS 


2.1  Fundamental  Equations 

2.1.1  General 


Hie  derivation  of  the  fundamental  equations  contained  herein  parallels  the 
derivation  of  the  steady- state  turbulent  boundary  layer  equations  for 
gaseous  systems.  The  effects  of  particles  are  included  with  the  assumption 
that  the  particles  and  gas  are  in  dynamic  and  thermal  equilibrium,  i.e.,  no 
velocity  or  temperature  lags.  Hie  particles,  therefore,  are  not  treated 
as  discrete  droplets,  but  as  chemical  species  which  have  properties  corre- 
sponding to  their  respective  liquid  or  solid  phases.  In  determining  the 
mass  transfer  of  the  particles  resulting  from  turbulent  mixing,  it  is 
assumed  that  the  particles  do  not  follow  the  random  motions  of  the  gas,  but 
due  to  their  inertia  diffuse  more  slowly  than  the  gas.  The  derivations 
that  follow  pertain  to  a two-dimensional  flow  model  and  later  are  expanded 
to  include  an  axisymmetric  flow  system. 


2.1.2  Global  Continuity  Equation 


Consider  a rectangular  control  volume  located  within  the  mixing  region  as 
shown  in  Figure  2.  It  is  assumed  that  variations  in  the  z-direction  are 
negligible.  The  flow  entering  the  elemental  volume  in  the  x-direction 
consists  of  a gas-particle  mixture  with  an  average  axial  velocity  u.  All 
the  flow  properties  are  assumed  to  be  time- averaged  values  consistent  with 
turbulent  theory  (8).  The  total  mass  flow  entering  the  volume  in  the 
axial  direction  is 


/OUAj  + tnpUAf, 
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vhere  /O  and  are  the  densities  of  the  gas  and  particles,  respectively 
and  Ay  is  the  flow  area  occupied  by  the  gas  and  that  by  the  particles 


Then 


fiif  * Af,  - Ay  AZ 


Assuming  Af>r* I,  and  letting  /o  equal  the  weight  of  particles  per  unit 


volume  of  gas,  yields 

(/%+/*,}  U A.y  A Z = /OtfAyAZ 

for  the  total  axial  mass  flow  into  the  volume  where 

/°  = /y  * ^ 

The  axial  mass  flow  leaving  the  volume  is  given  by 
jOt/AyAZ  + /^uAyAz)  AX 

The  net  change  in  axial  mass  flow  is 

-| -(/*u)  AXAyA? 

Similarly  in  the  y-direction  the  net  change  is 

2_  (f>zr)  AX  Ay  A Z 

Conservation  of  mass  yields 


■55  O")  * r 0 


for  the  steady- state  global  continuity  equation. 


2.1.3  Species  Continuity 


itions 


The  species  continuity  equations  are  derived  from  mass  balances  for  each  of 
the  chemical  species.  The  mass  balance  for  component  "i"  of  the  mixture 
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flowing  through  the  control  volume  consists  of  terms  involving  (a)  the  net 
mass  change  due  to  the  flow  entering  and  leaving  the  volume  and  (b)  production 
(or  depletion)  of  the  species  as  a result  of  chemical  reactions.  The  net 
transfer  of  component  i in  the  x-direction  due  to  the  fluid  motion  is 

where  is  the  mass  fraction  of  component  i (mass  per  unit  mass  of  mixture) 
and  is  its  mean  axial  velocity.  Similarly  in  the  y-direction  the  net 
transfer  is 

- OC/  ir.  ) 

The  net  rate  of  production  of  component  i is  given  by 


OUj  AX  Ay  Air 


A mass  balance  for  component  i yields 


According  to  theory  (9)  the  term  may  be  expressed  in  terms  of 

diffusion  and  convective  fluxes.  Thus 

p C u.  - hn,  + /o  c.u 


where  M ^ is  the  diffusion  flux  resulting  from  concentration  gradients 
and  P c u is  the  convective  flux  due  to  the  bulk  motion  of  the  fluid. 
By  bulk  motion  is  meant  the  average  macroscopic  motion  of  the  flow  in 
which  the  axial  bulk  velocity  u,  is  defined  by 


“ = Zw 
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vhere  h is  the  total  number  of  species  comprising  the  mixture.  For 
molecular  diffusion  the  flux  th4)(  , may  be  expressed  as 


' -*4  W 


(5) 


where  Di  is  the  molecular  diffusion  coefficient  for  component  i.  For  tur- 
bulent diffusion  it  has  been  shown  theoretically  (9)  that  a similar  expres- 
sion is  valid  with  the  molecular  diffusivity  replaced  by  its  turbulent 
counterpart . Thus , 

. C- 

- -S’  Tt 

(6) 

where  €pi  is  the  eddy  mass  diffusivity  of  component  i.  Similarly  in  the 
y-direction 


= *PJ 


3c; 


'/  ~PJ  3/  (7) 

Substitution  of  equations  3,  6,  and  7 into  equation  2 yields 

a (/°Cju)  2(/°CV)  _ ) f > (/*€„' V? ) ti,. 

dy  ~ dx  -*Y 

Expanding  the  above  equation  and  substituting  the  global  continuity 
equation  yields 


bx 


3 / 


UJj 


It  is  generally  assumed,  based  on  an  order  of  magnitude  analysis,  that  the 
first  term  on  the  right  side  of  the  above  equation  is  negligible.  Thus 
the  species  continuity  equation  for  component  i becomes 


(8) 
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She  net  rater,  of  production  of  the  species  are  determined  from  a considera- 
tion of  the  chemical  aspects  of  the  problem.  A detailed  discussion  of  the 
determination  of  li >.  is  presented  in  Section  4 herein. 

The  species  continuity  equation  developed  above  is  generally  assumed  valid 
for  gaseous  species,  with  the  eddy  mass  diffusivities  taken  to  be  identical 
for  all  species.  However,  for  particles  further  consideration  must  be  given 
to  the  turbulent  transport  mechanism.  The  turbulent  motions  within  a gas 
are  usually  conceived  as  macroscopic  "lumps  or  eddies"  of  fluid  moving 
randomly  from  one  region  to  another.  For  gas-particle  systems  it  seems 
unlikely  that  the  particles,  due  to  their  inertia,  are  capable  of  following 
the  random  fluctuations  of  the  gas.  Consequently,  the  net  transfer  of  part- 
icles from  one  region  to  another  is  likely  to  be  considerably  less  than  that 
of  the  gas.  This  difference  in  the  turbulent  transport  of  gas  and  particles 
may,  in  theory,  be  accomplished  by  employing  different  values  for  the  eddy 
mass  diffusivities  of  the  two  components.  Thus,  in  the  present  analysis  two 
values  of  the  eddy  mass  diffusivity  are  employed  - one  for  the  gas  and  the 
other  for  the  particles.  A detailed  discussion  of  the  numerical  values  to 
be  employed  for  the  eddy  mass  diffusivities  is  presented  in  Section  3 herein. 

2.1.4  Momentum  Equations 

The  axial  momentum  equation  is  derived  by  equating  the  net  change  in  axial 
momentum  of  the  fluid  flowing  through  the  control  volume  to  the  net  axial 
forces  acting  on  the  control  boundaries.  Thus 


*(/°«*) 


a (/*ua)  ?>p 


*v 


ax 


2>X 


(9) 
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where  7T„  is  the  shear  stress  which,  for  turbulent  flow,  may  be  expressed 
in  terms  of  the  velocity  gradient  and  the  eddy  momentum  diffusivity.  Thus 

Jxy  ' 3/ 

Expanding  equation  9 and  substituting  equation  1 yields 


(10] 

Employing  the  usual  boundary  layer  approximations  and  an  order  of  magnitude 
analysis,  the  y-direction  momentum  equation  reduces  to 


¥ --  ° 


2.1.5  Energy  Equation 

The  energy  equation  is  derived  following  the  methods  of  Schlichting  (10) 
and  Vasiliu  (3).  For  an  observer  moving  with  the  flow  the  First  Law  of 
Thermodynamics  may  be  stated  as  follows: 


[. 


rate  of  increase 
of  internal  energy 


C 


'work  done  by  the  system 
rate  of  heat  added!  + boundaries  on  the 

to  the  system  ; 


surroundings 


The  net  change  in  internal  energy  of  the  gas-particle  fluid  flowing  through 
the  control  volume  may  be  expressed  as 


dC  - (^| f + 


(11 


or 


dF  = /° 


Pe 

Pt 


AX  A/ A 2 


(12 
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where 


is  the  substantial  derivative  for  steady-state  and  e is  the  internal  energy. 

From  the  definition  of  enthalpy,  h*e  + r%,,  the  substantive  derivative  of 
the  internal  energy  may  be  written  as 

Pe  _ £h  __  P_(P/p)  (1 

pt  ~ pt  pt: 

where  the  enthalpy  of  the  mixture  is  given  by 

h=  fc.h;  - jj.  ( (Cr-dT  0 

The  heat  of  formation  i ih ^ , is  evaluated  at  the  reference  temperature, 
0°R.  The  substantive  derivative  of  enthalpy  is 


Ph  _ 7 (r  y. 

~DT  ~ IS  * Pt  ^ ' Dt 


Substitution  of  equations  14  and  1 6 into  12  yields 

- U ^)(C  &&  + h .££*')  - /O  AXAyA 2 

^ Zi  " Dt  M pt  1 pt 
L /-) 

or 

d£  L \Cj g £T  + K0L)  - /» 

, >=(  J 


dhx  _ 

Jr  “ ^ 


= C„ 


hc° 


Thus,  the  net  change  of  the  gas-particle  internal  energy  becomes 
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nr 

pt 


p (%>) 

r pt 


The  heat  transferred  or  added  to  the  control  volume  is  assumed  to  te  that 
resulting  from  the  following  processes:  (a)  conduction  and  diffusion,  (b) 

frictional  heating,  and  (c)  chemical  reactions.  The  net  heat  transfer  in 
the  y-direction  due  to  conduction  and  diffusion  may  be  expressed  as 


where  eH  is  the  eddy  diffusivity  for  heat  and  v?  is  the  y-direction 
diffusion  velocity  of  component  i defined  by 


xr:  - - 


*Y 


Thus  the  net  heat  transfer  term  becomes 


jQy  z (PCP**Xy) 

Expanding  the  second  term  in 
continuity  equation  yields 


the  above 


equation  and  substituting  the  species 


d<?,  = 47)  - "*) hj 

^ 7*7 

AS! 


The  frictional  heating  term  is  assumed  identical  to  that  generally  employed 
in  boundary  layer  analyses.  Thus 


aqf  --  fit,  (fy) 

(19) 
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Th.  heat  release  due  to  chemical  reaction  is  accounted  for  by  Including  the 
heat  of  formation  In  the  enthalpy  term  for  each  component . Thus  the  enthalpy 

of  the  mixture  is 


= 2 c-  f K 6T ' * 

>*  i ^ 0 


(20) 


•me  net  rate  of  vork  done  by  the  control  volume  surface  forces  on  the 
surroundings  Is  shovn  by  Schlichtlng  to  be  given  by 


(21) 


Employing  equations  17,  18,  19  and  21  the  energy  equation  becomes 


+ />*(! ff  ' t*'*  * V V 

Jtl  -**' 


The  second  term  on  the  right  side  of  the  above  equation  may  be  simplified 
utilizing  the  global  continuity  equation.  Thus 


P 
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Thus  the  energy  equations  becomes 
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2.1.6  Equation  of  State 

The  equation  of  state  for  the  gas-particle  mixture  is  developed  based  cn 
the  assumptions  that  (a)  the  gas  obeys  the  perfect  gas  law  and  (b)  the 
volume  occupied  by  the  particles  is  negligible.  Consider  a unit  mass  of 
mixture  comprising  a mass  Ciy  , of  gaseous  products,  and  a mass  CUp  of 
particles.  The  total  volume  of  the  mixture  is  given  by 

f°3  P 

where  and  are  the  densities  of  the  particles  and  gas,  respectively. 
Assuming  that  the  partial  volume  occupied  by  the  particles  is  negligible 
the  expression  for  the  total  volume  becomes 


Hie  equation  of  state  for  the  gaseous  phase  is 


P 


where  Rg  is  the  gas  constant  for  the  gaseous  components,  which  may  be 
expressed  in  terms  of  the  universal  gas  constant  and  the  gas  molecular  weight. 


*5  = 


where 
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The  is  perked  only  Pc,  the  6—  specie,  - * belts  the  — 
number  OP  saseous  components.  *»  total  density  - «-  — is  obtained 
from  equations  23-26.  Tiius 






r l'1 

p 

p ' KT 

[ 1 ] 

2.!.7  n+.emate  Format  of  the  Fundamental  Equations 

^ governs  equations  may  be  expressed  in  alternate  Ponss  employing 
dimepsioniess  ratios  involving  the  eddy  diPPusivlties.  » dePXPitiop  the 
turbulent  Leuis  and  Prapdtl  papers  axe 


/ 

Lej  c 2r— 


pb  = -£*- 


. (23) 
(29) 


vhere  „ and  C,  are  the  density  and  speciPic  heat  oP  the  gas-paxticXe 
Blrture,  respectively.  Employing  these  dimensionless  numbers  the  Pu 
equations,  in  axisymmetric  coordinates,  become 


Global  Continuity: 

5 (Pur') 

^ *>- 


Species  Continuity. 


/*"&  * ^TF 


3 I£i)  t UJJ 

TT^  Pb  T>r  ' 


(31) 


Momentum  Equation: 


P“  If  + ^ 


Du  _ + _L  l~(/°*vlr/’  T?) 

3^  ‘ ax  t*  *h 
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Energy  Equation: 

(^u  If  = "5?  r 'JpT^Cr'ri0p?^ 

Equation  of  State:  <'*1 


/*  = 


[ ! V*  I 

Ty  ' L >*=/  J 


The  exponent  for  the  terms  h* , makes  the  equations  applicable  for  either 
two-dimensional  (S=o)  or  axisymmetric  flow  (Sri).  A system  of  (4  + « ) 
simultaneous  equations  are  therefore  available  for  determining  the  unknown 
parameters  u,  v,  r,  c ■ , *&„■,  and  . Thus,  a solution  may  be  obtained 

if  appropriate  expressions  for  the  parameters  and  are  available 

in  terms  of  the  remaining  unknown  parameters.  Such  expressions  will  be 
developed  later  herein;  however,  for  the  present  these  parameters  will 
be  retained  in  their  present  form  with  the  assumption  that  suitable 
expressions  are  available. 


Transformation  of  the  Fundamental  Equations 


In  obtaining  a solution  it  is  convenient  to  transform  the  above  equations 
utilizing  the  von  Mises  transformation.  Instead  of  the  exisymmetric 
coordinates  x and  r (or  the  cartesian  coordinates  x and  y)  the  stream 
function  V , is  introduced  as  one  of  the  independent  variables  replacing 
the  r- coordinate.  The  x- coordinate  is  not  affected  by  the  transformation. 
By  definition  the  stream  function  is 


m. 


= /=>  “ r 


mr  ■ 


= Z^A 
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which  when  substituted  into  the  global  continuity  equation  (equation  30 ) 
yields 


i_  ( 2Jt\  _ -5-  l^-\  = 

3X  \ 3 h ) 'bfr  \ ! 


O 


Thus,  the  continuity  equation  is  satisfieu  identically  since 


>**>■  ' ~3!r>X 


The  transformation  of  the  remaining  equations  is  accomplished  by  noting, 
from  the  chain  rule  of  the  calculus,  that 

(3t\  _ fi£\  (M\ 

(*►■/*  " Jjr 

“ ft),  ’ m,  m,  • m. 

- (ii  l + (£  l 

where  f represents  any  of  the  dependent  variables.  Hence,  the  transformed 


equations  become 


Species  Continuity: 


2C;  _ 5 ( <£>  Les  jo3 u h**  ^C-A ) y.  J±!A 

~3X  ~ ’3?'  Pp  \ 


(37) 


Momentum  Equation: 


3U  _ _ / eip. 

2K  ~ .ou  dX 


+ 


i_r 


~s!AL  ) 
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(38) 
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The  finite  difference  technique  of  solving  differential  equations  consists 
of  replacing  the  partial  derivatives  by  finite  difference  ratios,  which  is 
equivalent  to  physically  replacing  the  continuous  flow  field  by  a network 
of  finite  elements.  Consider  the  flow  field  illustrated  in  Figure  1+  and 
assume  it  to  be  divided  into  a number  of  elements  having  grid  sizes  of  AX 
and  AtfJ  . If,  in  addition,  it  is  assumed  that  values  of  all  the  flow 
properties  are  known  along  the  vertical  line  designated  m,  referred  to  as 
"front  m”,  it  is  then  possible  to  calculate  the  flow  properties  at  front 
m + 1 in  the  following  manner. 


D2-36251-1 


Energy  Equation: 


37~  _ ! dP 

~ pCpd*  Co 


' 4 = 1 i Ti 

In  order  to  obtain  the  solution  in  the  physical  plane  (x,  r),  it  is 
necessary  to  transform  the  ^-coordinate  to  the  r-coordinate.  From  the 
definition  of  the  stream  function,  equation  35>  the  r-value  corresponding 


to  any  ip  -value  is  given  by 


+ (<?+i) 


where  r is  the  r-value  corresponding  to  the  reference  streamline,  ^ = 0. 

Shown  in  Figure  3 is  the  orientation  of  the  x-r  coordinate  system  and  the 
location  of  streamlines  for  a typical  problem. 

2 . 3 Development  of  the  Finite- Difference  Equations 
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Figure  4 - Finite  Difference  Grid  Network 


l.  i i , _ 


7or  the  grid  point  denoted(m  +1,  n)a  Taylor  series  expansion  of  the 
function  f in  terms  of  x- derivatives  yields 


f - -f 


AX  (¥)  + 

\*x 


AX*[ >V\ 

*!  \*x*L,t 


+ 


where  f may  be  taken  to  be  any  dependent  variable,  e.g.,  u,  t,  e,  etc. 
If  AX  is  taken  sufficiently  small  the  terms  of  higher  order  than  unity 
may  be  neglected,  thus  the  x-derivative  becomes: 


( 


If) 

2XJ, 


'h*,  h 


'hi,* 


AX 


The  above  equation  is  termed  the  forward-difference  approximation  for  the 
x-derivative.  In  a similar  manner  backward  - difference  and  average  - 
difference  approximations  may  be  derived  (7)*  The  resulting  backward- 
and  average-difference  approximations  are,  respectively 


/! 


vxj, 


7At<  H , 


fyj,n 


AX 


(■ 


'If) 

zxl 


tn.h 


"hn+l , H , 

^/]X 


Following  the  analysis  of  Wu  (7)  the  forward-difference  approximations 
are  used  herein  for  the  x- derivatives,  and  the  average- difference  approxi- 
mation for  the  (^-derivatives.  The  selection  of  these  approximations  is 
based  on  physical  and  stability  considerations  as  discussed  by  Wu.  In 
the  present  analysis,  the  following  approximations  are  employed  for  the 
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Substituting  the  above  expressions  for  the  derivatives  appearing  in 
equations  37-39  and  rearranging  terms,  yields 
Species  Continuity: 

C-  - C-  + (r.  - ZC.  + C.  ) (l 5) 

+ ^ V.1 

* ('oah%>  [ - (r^J,  * 

Momentum  Equation: 

u --  - (M3)  idf)  + AL(/>'«S%) ) 

*#"  ( />u  jtofUxL  a#  ' 

(46; 
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Energy  Equation: 


* 


T - T + I JP)  t-  — (.PU  tXF)  PrO°€v) 

Tf»+ 1,»~  I*."  ( JfiC-r  dK  L.n  t 9AV'*?r'-  '•».*»!.  A 

+ - *T»,*  + T^->)  +[('a*rh,—r  O0*rLs*$(%i**r  r* 


(47) 


&X  f/o'c/h1* *J\  \r  ia.(c.  - c.  \ fr  -r  ) — -4*-  Yftc/.  A.) 

WU'*"  W ,A  ^ 

As  can  be  seen,  the  above  equations  form  a set  of  algebraic  expressions  in 
which  the  dependent  variable  at  the  m + 1 front  is  expressed  explicitly 


in  terms  of  the  known  properties  at  front  m.  The  calculation  of  the  flow 
properties  throughout  the  grid  network  may  then  be  carried  out  by 
successively  applying  the  equations  to  each  front.  However,  it  is  first 


necessary  that  the  boundary  conditions  be  specified. 


2.4  Initial  and  Boundary  Conditions 

In  starting  the  calculations  for  a specific  problem,  it  is  necessary  that 
the  properties  along  the  initial  front  be  known  at  each  grid  point  - 
these  properties  are  termed  the  initial  conditions.  Hie  propeities  may 
be  constant  or  variable  dependent  on  the  problem  under  consideration.  For 
example,  it  may  be  uesired  to  assume  constant  properties  for  the  rocket 
exhaust  jet  except  in  the  immediate  neighborhood  of  the  nozzle  wall  where 
the  Jet  boundary  layer  exists.  Similarly,  the  flow  properties  in  the 
outer  stream  along  the  initial  front  may  be  variable  corresponding  to 
those  existing  within  a boundary  layer.  In  other  problems,  it  may  be 
desired  to  start  the  calculations  at  a point  downstream  where  the  mixing 
region  has  partially  developed  and  the  calculations  continued  for  further 
development  of  the  mixing  region.  For  any  of  these  situations,  it  is 
required  that  the  properties  be  known  at  each  grid  point  along  the  initial 
front  before  the  calculations  can  begin. 
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Once  the  initial  conditions  are  specified  the  flov  properties  at  each  grid 
point  along  the  second  front  may  be  calculated  except  at  the  grid  points 
corresponding  to  f - 0 and  Iff  - These  streamlines  are  the 

boundaries  of  the  grid  network  within  which  the  mixing  region  lies.  In 
order  to  obtain  a complete  solution  for  a given  flow  problem,  it  is 
necessary  that  the  flow  properties  along  the  boundaries  be  known.  These 
properties  are  termed  boundary  conditions.  It  is  the  usual  practice  in 
free  jet  problems  to  assume  the  flow  to  be  inviscid  in  the  areas  outside 
the  mixing  region.  The  boundary  conditions  for  such  problems  are  therefore 

specified  accordingly. 

For  problems  in  which  the  secondary  flow  is  contained  within  a duct, 
however,  three  regimes  of  flov  may  be  encountered  each  dictating  different 
boundary  conditions.  For  example,  consider  the  flow  field  shown  in 
Figure  5,  illustrating  the  different  flow  regimes.  The  first  regime 
begins  at  the  exit  plane  of  the  nozzle  and  ends  at  the  point  where  the 
inner  edge  of  the  mixing  zone  reaches  the  jet  centerline.  The  boundary 
conditions  for  this  regime  may  be  taken  to  correspond  to  those  for 
inviscid  flow.  The  flow  field,  therefore,  resembles  that  of  free  jet 
problems . 

In  the  second  regime,  the  mixing  region  widens  until  at  some  point  the 
outer  edge  of  the  mixing  zone  reaches  the  duct  wall,  whereupon  the  third 
regime  commences.  The  technique  for  determining  the  outer  boundary  condi- 
tion in  the  second  regime  is  identical  to  that  in  the  first,  i.e.,  the 
secondary  flow  outside  the  mixing  region  is  assumed  to  be  inviscid.  The 
inner  boundary  conditions  for  the  second  and  third  flow  regimes  may  be 
determined  from  the  physical  characteristics  of  the  flow.  Along  the  jet 
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Flow  Regimes  for  Confined  Mixing 


centerline,  the  flew  properties  change  as  a result  of  the  mixing  process. 
However,  from  symmetry  it  may  be  concluded  that  on  the  centerline,  the 
radial  gradients  of  all  properties  must  be  zero.  Ibis  condition  may 
therefore  be  employed  as  the  inner  boundary  condition  for  the  second 

regime . 

In  the  third  regime,  the  inviscid  outer  boundary  assumption  is  no  longer 
valid  since  the  mixing  region  extends  across  the  entire  duct.  However, 
the  outer  boundary  condition  for  this  regime  may  be  established  by 
assuming  that  the  heat  transfer  and  shear  stress  at  the  duct  wall  are 
negligible.  Thus,  from  the  rate  equations  for  shear  stress  and  heat  transfer 

it  is  seen  that 

/l")  - 0 

I ’bV'hfjan 

(%l,  * • 

In  addition,  since  mass  cannot  be  transferred  across  the  duct  wall,  we  have 


Thus,  the  inner  and  outer  boundary  conditions  for  the  third  regime  are 
identical.  Summarizing,  the  boundary  condixions  for  the  first  flew  regime 


may  be  expressed  in  the  following  form. 
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The  above  equations  are  applicable  for  determining  both  the  inner  and 
outer  boundary  conditions.  In  the  second  regime,  the  boundary  conditions 
for  the  inviscid  region  are  as  above;  whereas,  the  inner  boundary  con- 
ditions are  given  by 


($Cj)  r 

(2“)  = 

fir) 

1 ’tV'lyso 

[ wJr-p 

= O 


(51) 
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For  the  third  regime,  the  boundary  conditions  are  expressed  in  the  form 
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It  should  be  noted  that  in  the  above  discussion  it  was  assumed  that  the 
inner  edge  of  the  mixing  zone  reached  the  duct  centerline  before  the  outer 
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edge  reached  the  duct  wall.  In  certain  situations,  however,  these  conditions 
may  be  reversed,  the  outer  edge  reaching  the  duct  wall  prior  to  the  inner 
edge  reaching  the  centerline,  in  which  case  the  aforementioned  boundary 
conditions  are  altered  accordingly. 
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[ 0 TURBULENT  TRANSPORT  COEFFICIENTS 


3.1  General 


in  order  to  obtain  a solution  of  the  finite-difference  equations,  it  is 
necessary  that  expressions  for  the  eddy  dlffusivities  be  established  in 
terms  of  known  quantities.  Unfortunately  the  eddy  dlffusivities  are  not 
fluid  properties,  as  are  their  molecular  counterparts,  but  parameters 
dependent  upon  the  fluid  motion,  e.g.,  level  of  turbulence,  velocity 
gradient,  density  gradient,  etc.  Since  turbulent  flow  is  characterized 
by  random  fluctuations,  the  derivation  of  accurate  theoretical  expressions 
for  the  eddy  dlffusivities  requires  a statistical  analysis  such  as  that 
initiated  by  Taylor  (11).  However,  the  statistical  approach  at  the  present 
time  has  not  proven  practical  (12).  Consequently,  one  must  rely  on  the 
semi-empirical  approach  for  determining  expressions  relating  the  eddy 
dlffusivities  to  the  flow  properties. 

3.2  Incompressible  Theory 

Most  of  the  semi-empirical  theories  to  date  are  extensions  or  modifications 
of  Prandtl 1 s mixing  length  concept.  From  Prandtl's  initial  hypothesis,  it 
was  found  that  the  shear  stress  for  incompressible  turbulent  flow  can  be 

determined  from 


r= 


W 


vhere  9 is  Prandtl's  mixing  length  parameter,  assumed  to  be  independent 
of  the  y-coordinate  (l,  10).  Broussinesq,  had  earlier  postulated  that, 
in  analogy  with  the  shear  stress  lav  for  laminar  flov,  the  shear  stress  for 
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turbulent  flow  could  be  expressed  in  the  form: 


where  is  the  eddy  diffusivity  of  momentum,  often  referred  to  as 
eddy  viscosity,  being  analogous  to  molecular  kinematic  viscosity.  From 
equations  S3  and  5^  it  is  found  that 


p a ~*u 
A ~by 


In  an  attempt  to  obtain  more  realistic  values  of  €v  when  approaches 

zero,  Prandtl  modified  the  above  expression;  however,  the  resulting  expression 
is  difficult  to  use. 


A simpler  expression  for  eY  was  subsequently  developed  by  Prandtl  based 
on  the  experimental  data  of  Reichardt.  Assuming  that  e*.  was  constant 
across  the  mixing  zone,  i.e.,  independent  of  the  y-coordinate,  Prandtl 
postulated  that  for  free  jets 


tfy  — K b ( ~~  ) 


where  k is  an  empirical  constant,  b is  the  width  of  the  mixing  region,  and 
(umax"^min)  is  the  difference  between  the  maximum  and  minimum  velocities 
in  the  mixing  region.  Furthermore,  it  was  found  that 

b = C X 


so  that  Eqn.  5 6 becomes 
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where  the  product  kc,  is  replaced  by  another  proportionately  constant  K. 
From  Reichardt's  experiments  of  incompressible  free  jets  (Umin  **  0),  it 
was  determined  that  the  width  of  the  mixing  zone  varied  according  to  the 
relationship 

b - 0.0*18  X (5* 

The  resulting  expression  for  which  best  correlated  the  data  was 


€v  - O.  00/37  X UMa7( 


An  extension  of  the  above  expression  is  often  employed  for  the  case  of  two 
moving  streams,  thus 

~ 0-00137  X (U^ax-  U^ih)  (6o) 

An  alternate  form  of  Equation  57  is  frequently  used  in  correlating  experi- 
ment and  theory  in  terms  of  the  similarity  parameter,  thus 


= 


(Syria  X 


where  o~  the  similarity  parameter,  is  a constant  to  be  determined  from 
experimental  data.  Physically,  tr  is  a parameter  related  to  the 
spreading  rate  of  the  mixing  region  - larger  values  of  <r  corresponding 
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to  smaller  spreading  rates.  From  the  data  of  Relchardt,  employed  In 
obtaining  Equation  60,  the  corresponding  value  of  or  vas  determined  to  be 

13-5- 

still  another  expression  for  <F„  is  that  given  by  Pai  (12),  which  in  its 
most  general  form  is 

ev  - £0  ( c 1-  x/l)  (62) 

where  <fe  is  referred  to  as  Reichardt's  coefficient,  C is  a constant  which 
accounts  for  a virtual  origin  of  the  mixing  region,  L is  a reference  length, 
and  n is  a number  having  a value  between  0 and  1 depending  on  the  degree  of 
mixing.  The  constant  C may  be  assumed  to  account  for  the  bound  ry  layers 
that  build  up  in  the  two  streams  prior  to  the  point  of  initial  contact 
(x  = 0)  which  produce,  in  effect,  a mixing  region  of  finite  width  at  x = 0. 

For  the  case  where  C - 0 and  n • 1,  as  generally  is  assumed,  Equation  62 

becomes 

= e*  x/l  <63) 

Equating  Equations  6l  and  63  and  solving  for  £0  yields 

^ L ( U*ruh ) (64) 

= ^ 

From  extensive  data  of  free  Jets  exhausting  Into  a qulesclent  atmosphere, 
the  value  of  <r  for  Jet  Mach  numbers  less  than  unity  has  been  relatively 
veil  established  as  being  cr  * 12.  For  cases  vhere  the  receiving  medium  is 
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also  moving  the  value  of  O"  is  found  to  depend  on  the  velocity  and  density 
ratio  between  the  two  streams  as  discussed  below. 


3.3  Compressible  Theoi 


In  the  development  of  the  expressions  for  incompressible  flow  it  was 
assumed  that  €¥  was  constant  across  the  mixing  zone,  an  assumption  which  has 
been  substantiated  by  extensive  experimental  data.  However,  for  compressible 
flow  involving  large  temperature  and  velocity  differences  between  the  two 
streams,  it  has  been  postulated  that  is  dependent  not  only  on  the  x- 
coordinate,  but  the  y- coordinate  as  well  (13)*  Several  mathematical  models 
for  the  compressible  eddy  viscosity  have  been  proposed  by  various  investi- 
gators (lU,  15,  16),  which  for  the  most  part  are  more  complicated  and  con- 
sequently more  difficult  to  use  than  the  incompressible  models.  In  addition, 
since  only  limited  experimental  data  are  available,  none  of  the  proposed 
compressible  models  has  been  satisfactorily  verified.  Consequently  it  has 
been  the  general  practice  to  employ  the  incompressible  model  for  £v 
and  adjust  the  empirical  constants  to  account  for  the  effects  of  compressi- 
bility. 

The  most  frequently  used  parameter  in  correlating  experiment  and  theory  is 
the  similarity  parameter  cr  f defined  in  Equation  6l.  Experimentally 
determined  values  of  cr  for  free  Jets  exhausting  into  a quiescent  atmosphere 
have  been  obtained  for  jet  Mach  numbers  ranging  up  to  approximately  3*0  (15)* 
Figure  6 illustrates  the  variation  of  C r with  jet  Mach  number  for  a jet 
exhausting  into  still  air.  The  solid  line  faired  through  the  experimental 
data  was  originally  presented  in  Reference  15 . The  dashed  line  illustrates 
the  theoretical  variation  of  cr  with  Mach  number  as  given  by  Korst  (17)*  The 
majority  of  the  experimental  data  were  obtained  with  air  jets  having 
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temperatures  approximately  equal  to  that  of  the  ambieht  air. 


exception  is  the  data  of  Anderson  and  Johns,  which  were  obtained  with  a 
solid  propellant  rocket  exhausting  into  still  air.  *or  Mach  n 
than  unity,  the  value  of  ^ has  been  relatively  veil  established  as  being 
o-  = 12.0.  For  Mach  numbers  greater  than  unity,  the  data  illustrate  that 

increases  considerably. 


For  the  turbulent  mixing  between  two  compressible  coaxial  streams,  analytical 

attempts  (18)  have  resulted  in  relating  the  tvo-stream  CT  -value 

velocity  ratio  between  the  two  streams.  *>e  equivalent  single  stream  value 

is  obtained  fromKorst's  empircal  relationship 


a-j.  - it  + *•  I*8 


m 


simile  stream  Mach  number  rented  to  the  Jet 
where  Mia  is  an  equivalent  single 

stream  Crocco  number  as  shown  in  Fig.  7.  Also  ahovn  in  Fig.  7 1=  the  ratio 
of  the  aforementioned  <T  -values  as  a function  of  the  velocity  ratio  bet 
the  two  streams;  indicating  that  the  two  stream  .-value  increases  with 

TMs  irend  has  also  been  established  from  expert- 
increasing  velocity  ratio.  This  trena  n 

mental  data  (12). 
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Since  experimental  data  for  predicting  two  stream  cr  -values  are  limited, 
the  above  technique  for  evaluating  cr  has  been  adopted  herein  with  the 
exception  that  the  equivalent  single  stream  a,  value  is  obtained  from 
Fig.  6 corresponding  to  the  equivalent  single  stream  Mach  number, 
corresponding  value  of  Reichard’s  coefficient  is  then  detained  from  Equation 


Velocity  Ratio  - U9/Up 

Theoretical  Effect  of  Velocity  Patio  on  Similarity  Parameter  Rrtlo 


It  should  be  noted  that  in  utilizing  any  of  the  above  €v  expressions, 
difficulties  may  arise  resulting  from  the  method  employed  in  correlating 
experiment  and  theory.  Hie  accepted  method  of  establishing  the  validity  of 
a particular  model  is  to  employ  a specific  theoietical  analysis  and  several 
trial  values  of  the  empirical  constant  in  the  expression,  until  agreement 
is  obtained  between  the  measured  and  calculated  velocity,  temperature,  and 
concentration  profiles.  Unfortunately,  the  resulting  €v  expressions 
reflect  any  peculiarities  of  the  particular  theoretical  analysis  used  in 
the  correlation,  and  e'jnsequfntly  should  be  u with  renf-rmticn  for  ot£*>r 
analyses . 


3.4  Turbulent  Lewis  and  Prandtl  Numbers 


The  turbulent  Lewis  and  Prandtl  numbers  are  the  turbulent  counterparts  of  the 
corresponding  numbers  based  on  molecular  properties.  By  definition 


u.  - f* 

X € H 


Pb  - 


Another  dimensionless  ratio  frequently  employed  is  the  Schmidt  number  defined 


<7  - 

5ci  ' Leu 


where  the  subscript  i infers  that  the  various  species  may  have  different 
eddy  mass  diffusivities. . In  early  mixing  analyses  it  was  often  assumed, 


in  analogy  with  the  molecular  ratios,  that  Pr  = Sc  = Le  = 1;  an  assumption 


which  Simplifies  the  governing  equations  considerably.  However,  recent 
experimental  data  indicate  that  these  parameters  may  differ  markedly  from 
unity  (l4,  19).  For  example,  from  the  experimental  results  of  Forstall  and 
Shapiro  (19),  it  was  found  that  the  turbulent  Pr  and  Sc  numbers  were  equal 
and  constant  throughout  the  mixing  region,  with  a value  of  approximately 
0.7.  From  the  experimental  results  of  Zakkay,  et  al,  (l4)  it  was  concluded 
that  the  turbulent  Schmidt  and  Lewis  numbers  ranged  from  0.3  to  2.3  and  0.4 
to  1.0,  respectively. 

In  the  present  analysis,  values  of  the  turublent  Prandtl  and  Lewis  numbers 
may  be  employed  other  than  unity,  however  they  remain  as  constants  independent 
of  the  x-y  coordinates.  For  gaseous  components  the  values  of  Pr  and  Sc 
most  frequently  cited  in  the  literature  vary  between  0.5  and  1.2.  It  was 
shown  in  Ref  . 20  that  variations  of  Sc  have  negligible  effects  on  the  radial 
profiles,  so  that  any  value  within  the  above  range  may  be  used  with  reasonable 

accuracy. 

-Hie  Lewis  and  Prandtl  numbers  for  the  solid  and/or  liquid  particles  in  gas- 
particle  mixtures  are  at  the  present  time  unknown.  The  experimental  data 
of  longwell  and  Weiss  (21 ),  however,  indicate  that  the  ratio  of  gas  to 
liquid  eddy  mass  diffasivities  is  approximately  2.  In  their  studies  liquid 
diesel  fuel  was  injected  into  a turbulent  air  stream  and  concentrations 
measured  at  various  axial  locations.  The  results  were  compared  to  similar 
data  obtained  with  naphtha  as  the  injected  fuel.  As  noted  by  the  authors 
the  naphtha  evaporated  readily  and  was  assumed  to  be  mainly  in  the  gaseous 
phase.  The  results  showed  that  the  ratio  of  the  gas  to  liquid  eD  varied 
from  approximately  1.2  to  2.0  for  gas  stream  velocities  ranging  from  200  to 

500  fps,  respectively. 
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4.0  CHEMICAL  MODEL 


4.1  General 

In  considering  the  chemical  aspects  of  problems  involving  the  mixing  between 
reactive  components,  two  choices  of  chemical  behavior  are  available j namely, 
equilibrium  or  nonequilibrium  chemistry.  Because  of  the  high  velocities 
encountered  in  rocket  exhaust  Jets,  it  appears  that  in  such  problems  non- 
equilibrium  chemistry  would  prevail  throughout  most  of  the  flow  field.  The 
most  accurate  analysis  should,  therefore,  be  based  on  nonequilibrium  flow 
including  the  complicated  chemical  equations  describing  the  reaction 
mechanism.  However,  in  view  of  the  complexity  of  such  an  analysis  and  the  re- 
sulting lengthy  calculations,  it  is  advantageous  to  use  equilibrium  chemistry 
whenever  that  choice  leads  to  accurate  results. 

In  the  discussion  that  follows  the  equations  for  both  equilibrium  and  non- 
equilibrium  chemistry  are  developed  for  the  hydrogen- oxygen  system  and  the 
results  compared.  From  the  comparison  of  temperature  profiles,  it  was  ascer- 
tained that  the  assumption  of  equilibrium  chemistry  was  valid  for  hydrogen- 
oxygen  systems.  Ibis  assumption  was  therefore  employed  for  such  systems 
and  extended  to  other  systems  as  discussed  below. 

4.2  Nonequilibrium  Chemistry 

For  nonequilibrium  chemistry  the  net  rate  of  production  of  the  species  vk,  is 
determined  using  the  methods  of  theoretical  reaction  kinetics.  For  example, 
consider  the  chemical  reaction  between  hydrogen  and  oxygen  forming  water 
vapor.  It  is  well  known  that  the  reaction  mechanism  between  these  molecules 
actually  consists  of  several  chain-type  reactions  including  the  following: 
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Hg  + OH  :=■  H20  + H 


H + Og  = OH  + 0 
Hg  + 0 ==■  OH  + H 

H+H  + M = Hg  + M 
H+OH  + M^=HgO  + M 
0 + 0 + M=Og  + M 


(TO) 


where  M is  a third  body  usually  taken  to  represent  all  the  molecules  in  the 
mixture. 


In  order  to  determine  the  net  rates  of  production  of  the  various  species 
entering  the  above  reactions,  it  is  necessary  to  apply  the  reaction  kinetics 
equations  to  each  of  the  reaction  equations.  In  order  to  simplify  the  present 
analysis,  however,  it  is  assumed  that  the  above  complex  chain  reactions  for 
the  Hg-Og  system  may  be  replaced  by  the  following  one-step  reaction  equation 

2H2  + 02  2H20  (71) 

where  kf  is  an  "overall"  reaction  rate  constant. 


Following  Penner  (22)  it  can  be  shown  that  the  net  rate  of  production  of  HgO, 
in  accordance  with  reaction  equation  71,  may  be  expressed  as 


u> . - z 


w,  K/°3  (— ] 

‘ l Wxl 


(£l 

K'pf?r/Tc3 

(72) 


where  k^.  is  the  forward  reaction  rate  constant,  Kp  is  the  equilibrium  constant 
and  W is  the  molecular  weight.  The  subscripts  1,  2,  3 correspond  to  the 
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species  HgO,  Hg  and  02,  respectively. 

From  further  application  of  the  reaction  rate  equation**  it  may  be  shown  that 
the  net  rate  of  production  of  Hj,  is  related  t0  that  °f  Hs°‘  ^ 

U,,  = ±4  (73 

From  conservation  of  mass  the  net  rate  of  production  of  u2  becomes 


= - ft  ^ 


(7*0 


•me  expreseion  for  the  reaction  rate  constant  kf  is  assumed  to  be  of  the  form 


5 7" 


(75) 


where  B is  the  frequency  factor,  E is  the  activation  energy,  and  Ms  a 
constant . 

Equations  (72)  through  (7k)  may  be  employed  in  solving  the  finite  difference 
energy  and  species  conservation  equations  provided  the  empirical  constants 
in  equation  (75)  may  be  established.  The  primary  difficulty  in  all  non- 
equilibrium problems  is  that  of  establishing  the  reaction  mechanism  and 
associated  expressions  for  the  reaction  rate  constants.  For  the  Sj-Oj 
system  the  complex  chain  reaction  equations  and  rate  constants  are  relatively 
veil  established.  It  is  possible  therefore,  to  establish  an  expression  for 
an  "overall"  reaction  rate  constant  in  the  following  manner.  Utilising  the 
chain  reactions  and  the  associated  reaction  rate  constants,  the  temperature 
distribution  throughout  a stream  tube  may  be  calculated  for  a range  of  initial 
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conditions  of  the  flow  parameters  employing  the  one-dimensional  analysis  of 
Ref.  23.  Then,  using  the  same  analytical  technique  and  the  one-step  chemi- 
cal reaction  a value  of  the  overall  reaction  rate  constant  may  be  deter- 
mined which  best  reproduces  the  previously  determined  temperature  distribution. 

The  above  technique  was  employed  for  a range  of  initial  conditions  correspond- 
ing to  those  anticipated  for  the  problem  being  considered  herein.  Typical 
results  are  shown  in  Fig.  8 wherein  are  presented  the  calculated  temperature 
distributions  resulting  from  the  use  of  chain  reactions  (the  solid  line) 
and  the  one-step  reaction  using  various  expressions  for  the  overall  reaction 
rate  constant  (the  broken  lines) . The  reaction  equations  employed  in  calcu- 
lating the  solid  line  were  those  presented  in  equation  70  herein.  The  value 
of  the  overall  reaction  rate  constant  was  varied  by  altering  the  frequency 
factor;  the  activation  energy,  obtained  from  Ref.  3,  was  assumed  constant 
at  1.6  x 10A  cal/gm-mole.  The  average  value  of  the  frequency  factor  which 
best  correlated  the  results  obtained  from  the  complex  reactions  was 
B - 1017.  Thus,  the  expression  employed  for  the  overall  reaction  rate 


constant  for  the  reaction  was  determined  to  be 


n n o.  a" 

= 16  T exr 


^ -i-H5 * to4  j 


(76) 


4.3  Equilibrium  Chemistry 

Employing  the  assumption  of  equilibrium  chemistry  the  net  rate  of  production 


The  w ± terms  appearing  in  the  finite  difference  equations  are  then  determined 
from  equations  72-74. 
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of  each  species  nay  be  determined  from  the  equation 


oj.  ~ 
A 


(<L,  - cu) 


where  and  are  the  concentrations  of  species  K inmediately  before  and 
after  the  chemical  reaction,  respectively.  The  problem,  therefore,  reduces 
to  that  of  calculating  the  equilibrium  composition  CK  , for  initial  concen- 
trations Ck  at  a given  pressure.  In  the  discussion  that  follows  the  equili- 
brium equations  are  developed  not  only  for  the  hydrogen-oxygen  system,  but 
in  addition,  for  systems  comprising  the  following  species:  HgO,  Hg,  0^  CC^, 
CO,  N2,  HC1,  Al,  and  AlgO^.  The  chemical  reaction  equations  are  assumed  to 


be  as  follows: 


2H2  + 02  2H20  (To; 

2 CO  + O2  =5=^=  2C0g  (79) 

2A1  + 3/2  02  = A1203  (6°) 

Consider  now  the  elemental  control  volume  shown  in  Fig.  2.  As  the  fluid 
traverses  the  distance  Ax,  the  concentrations  change  as  a result  of  the 
mixing  process.  If  it  is  assumed  that  the  fluid  entering  the  volume  is  in 
chemical  equilibrium,  then  at  x + A x the  fluid  is  not  in  equilibrium.  Since 
the  assumption  of  equilibrium  flow  implies  that  the  reactions  occur  at  an 
infinite  rate,  the  composition  at  x + A x may  be  determined  as  follows. 

From  an  atom  balance  the  following  relations  may  be  obtained  which  are 
valid  for  both  reactants  and  products: 


C - c -h  c 
S*  1 ' w , 


'll  . 
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c6  = ^6 

C7  - £7 


C„=  ^ 


i-  c ^ 

, c* 

* 3 

IV, 

> 2W, 

t 

s zws 

z 

*"> 

II 

* cs 

w* 

- ■ 

C7=  C, 


c, 


where  is  a pseudo-mass  fraction  defined  by  the  above  equations  and 
the  subscript  i = 1,  2,  - - -,  10  refers  to  the  species  HpO,  H~,  0 2,  CC^, 
CO,  Ng,  HC1,  any  inert  component,  Al,  AlgO^,  respectively.  Equations  8I-87 
are  also  valid  for  the  products  - the  cj  on  the  right  side  replaced  by  c-^. 

The  equilibrium  constants  for  the  reaction  equations  78-80,  expressed  in 
mass  fractions,  are 


4 = 


/ C W,  \* 

V w,cj 

(88) 

( <k 

(89) 

W \Vx  ~7/x 

7?  (rw„) 

(90) 
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where  the  fugacity  of  AI2O3  is  assumed  to  be  unity  (24),  and  Wm  is  the 
average  molecular  weight  of  the  gas  mixture.  Values  of  the  equilibrium 
constants  are  obtained  from  tables  of  Kp  at  various  temperatures  (25). 
Equations  81-90  are  seen  to  form  a system  cf  ten  equations  for  the  ten 
unknown  concentrations  c*.  The  solution  for  the  concentrations  comprises 
an  iteration  technique  discussed  in  Ref.  6.  The  net  rate  of  production 
terms  are  then  calculated  from  Equation  77  • 

4.4  Comparison  of  Equilibrium  and  Nonequilibrium  Chemistry 
In  order  to  illustrate  the  difference  between  equilibrium  and  nonequilibrium 
chemistry,  solutions  were  obtained  for  a typical  problem  employing  both 
chemical  models  for  the  hydrogen- oxygen  system.  Shown  in  Fig.  9 is  the 
temperature  profile  throughout  the  mixed  flow  field  at  an  axial  distance 
of  four  feet  downstream  from  the  exhaust  nozzle  of  0 LC^/LHp  rocket  motor . 
The  conditions  of  the  rocket  exhaust  and  secondary  air  stream  at  the  exit 
plane  of  the  nozzle  are  noted  in  the  figure.  As  can  be  seen,  the  difference 
between  the  two  cases  is  not  great,  thus  demonstrating  that  the  selection 
of  equilibrium  chemistry  provides  sufficiently  accurate  results  for  the 
given  conditions.  As  a point  of  interest,  it  is  worthwhile  to  note  that 
the  machine  computation  times  for  equilibrium  and  nonequilibrium  flow  were 
approximately  2 and  110  minutes,  respectively.  In  view  of  the  long  computa- 
tion times  required  for  nonequilibrium  flow,  and  since  the  assumption  of 
equilibrium  flow  yields  reasonably  accurate  results,  the  latter  chemical 
model  was  assumed  for  the  problem  under  study  herein. 
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5.0  COMPUTER  PROGRAM 


5.1  Capability 

Although  the  computer  program  described  herein  was  developed  specifically 
for  the  air-augmentation  problem  it  may,  however,  be  employed  for  other  re- 
lated problems.  For  example,  the  program  may  be  utilized  in  calculating 
the  mixing  region  for  free  Jet  problems,  where  the  ambient  atmosphere  is 
either  m-ving  or  at  rest.  In  addition,  calculations  may  be  performed  for 
problems  involving  the  mixing  and  combustion  of  fuel  injected  into  an  air 
stream,  with  application  to  ramjet  cctnbustor  studies.  Problems  involving 
chemically  reactive  species  may  be  solved  employing  either  equilibrium  or 
frozen  chemistry.  It  is  noted  that  the  program  is  limited  to  the  calcu- 
lation of  turbulent  flow  fields  - laminar  flow  problems  being  excluded. 

In  the  following  sections  will  be  described  the  input  options  available  to 
the  user,  the  important  features  of  the  program,  and  the  procedure  to  be 
followed  in  establisldng  the  input  data  for  typical  problems.  For  con- 
venience in  the  discussion  that  follows  a list  of  the  pertinent  computer 
program  notation  is  presented  in  Appendix  7«1. 

5.2  Pressure  or  Duct  Profile  Option 

In  air-augmentation  problems  it  is  generally  desirable  to  specify  the  duct 
geometry;  whereas,  in  free  Jet  problems  it  is  more  convenient  to  specify 
the  pressure  distribution.  In  the  computer  program  the  option  is  available 
to  specify  either  the  duct  profile  in  the  form  rD  = f(x),  or  the  pressure 
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distribution  in  the  form  p = g(x).  The  program  control  for  this  option  is 
termed  INDUCT  - having  values  of  0 or  1 for  pressure  distribution  or  duct 
geometry  input,  respectively.  The  input  of  p » g(x),  or  rD  = f(x),  is  in 
the  form  of  a table.  The  x-values  are  termed  XT AS  and  the  p-  or  r^-  values 
termed  PTAB  in  program  notation.  It  should  be  noted  that,  from  the  view- 
point of  machine  time,  it  is  advantageous  to  input  the  pressure  distribution 

whenever  possible. 

If  the  pressure  input  option  is  chosen  the  duct  radius  is  calculated  from 
the  initial  specified  duct  radius,  denoted  RDZERO,  and  conservation  of 
mass  of  the  flow  through  the  duct.  Conversely,  if  the  duct  input  option 
is  chosen,  the  pressure  at  each  front  is  calculated  from  thr-  initial  speci- 
fied pressure  PZERO,  and  conservation  of  mass.  Problems  involving  flee  jets 
are  treated  by  specifying  RDZERO  very  large  compared  to  the  nozzle  exit 

radius  RE. 

5.3  Reference  Streamline 

In  performing  the  finite-difference  calculations  it  is  necessary  that  the 
coordinates  of  a reference  streamline  be  specified,  i.e.,  = f(x>r)» 

Thus  the  inverse  transformation,  given  by  Eqn.  40,  may  be  performed.  In 
the  present  program  the  reference  qr -line  may  be  specified  as  coincident 
with  either  the  rocket  exhaust  plume  boundary  or  the  centerline  of  the 
flow  field  as  shown  in  Fig.  10.  In  the  former  case  the  coordinates  of  the 
plume  slipline  may  be  determined  utilizing  the  method  of  characteristics 
program  of  Ref.  5.  The  slipline  coordinates  are  input  to  the  program  in 
tabular  form,  denoted  as  XTAB  and  RTAB.  The  input  value  of  RZER01  is 
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accordingly  specified  as  being  slightly  larger  than  zero,  e.g.,  RZER01  «* 
0.001  is  usually  adequate.  One  problem  which  is  generally  encountered 
using  the  slipline  as  f ref,  is  as  follows.  As  the  calculations  proceed 
in  the  x-direction  the  mixing  zone  widens,  until  at  a certain  distance  the 
mixing  zone  crosses  the  ^>=  0 streamline.  Once  thiB  condition  occurs  the 
reference  y'-line  is  shifted  to  the  jet  centerline  and  the  calculations, 
continued,  or  if  desired  the  program  is  stopped.  The  control  parameter 
for  this  option  is  termed  ISET  - having  values  of  0 or  1 for  stopping  or 
continuing  the  program,  respectively.  For  cases  where  it  iB  desired  to 
specify  the  centerline  as  the  reference  ^ -line  the  value  of  RZER01  and 
RTAB  are  set  equal  to  zero  and  the  radius  of  the  jet,  respectively. 


5.4  Mixing  Zone  Definition 


The  width  of  the  mixing  zone,  as  defined  herein,  is  determined  from  the 
species  concentrations.  At  a given  "front",  the  calculation  of  all  flow 
properties  begins  in  the  mixing  region  and  proceeds  inward  and  outward 
until  the  edges  of  the  mixing  region,  as  established  from  the  following  cri- 
teria, are  determined.  At  each  grid  point  the  difference  is  calculated  be- 
tween the  concentration  of  species  i at  the  grid  point  in  question  and  the 


outer  boundary  condition  of  C^,  denoted  as  £luan^y  *-s  ca^~ 


culated  for  each  of  the  species,  i - 1-10,  and  the  maximum  difference 
determined.  The  maximum  difference  is  then  compared  with  an  input  quantity 
termed  TESMAX.  The  outer  edge  of  the  mixing  region  is  established  when 


(cjn  ' ^Wr)Mfx 


res  max 


(91) 


BgjEi/SfG 

NO  D2-36251-1 

REV  LTR 

sh.  63 

U3  4288-2000  REV.  6/64 


Similarly,  the  inner  edge  of  the  mixing  zone  is  established  when 


(cs„  - (92) 

where  CiJMIN  is  the  inner  boundary  condition  of  C^,  and  TESMIN  is  an  input 
quantity.  The  values  of  TESMIN  and  TESMAX  should  be  sufficiently  small  so 
that  the  mixing  zone  boundaries  determined  from  species  concentrations  are 
larger  than  the  corresponding  boundaries  determined  from  velocity  or  tem- 
perature profiles.  A value  of  0.00001  for  both  TESMIN  and  TESMAX  is  suf- 
ficient for  most  problems  - smaller  values  yield  wider  mixing  zones. 

5.5  Stability  Criteria 

The  major  difficulties  in  living  differential  equations  utilizing  the 
finite-difference  technique  involve  the  problems  of  stability  and  con- 
vergence of  the  solution.  A stable  solution  is  defined  as  one  where  small 
errors,  which  are  inherent  in  the  calculations,  do  not  increase  in  size. 

In  unstable  solutions  the  errors  increase  in  size  during  the  calculations 
and  eventually  become  so  large  that  the  results  are  noticeably  inaccurate. 

A solution  is  defined  as  convergent  if  by  decreasing  the  grid  size  ( dx  and 
A V)  the  calculated  velocity,  temperature,  and  concentration  profiles  are 
not  significantly  altered.  It  has  been  established  that  if  a solution  is 
stable  it  is  also  convergent  (7,  26). 

Buploying  Kar plus'  theory  of  stability  it  was  determined  that  for  the  pre- 
sent problem  the  species  continuity  equation  dictated  the  most  stringent 
stability  critsria  (5).  Thus  in  order  to  insure  a stable  solution,  and 
consequently  a convergent  solution,  the  stability  criteria  dictates 
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that 


W * (?/o*ur  evM-pt  ) 


where 

^ ( C f x/l  f (94) 

In  application,  the  right  side  of  Eqn.  93  is  evaluated  at  each  grid  point 
in  the  mixing  region  along  a given  front,  and  the  maximum  value  determined. 
For  stability  the  value  of  &ty  must  be  greater  than  or  equal  to  the  estab- 
lished nmTHimirn  value.  As  can  be  seen  the  value  of  the  right  side  of  Eqn.  93 
is  related  to  the  axial  distance.  Thus,  it  becomes  evident  that  if  constant 
values  of  Aty  and  4x  are  selected,  special  consideration  must  be  given  to 
the  selection  of  their  values  so  as  to  insure  stability  throughout  the  cal- 
culations. 

In  order  to  remove  the  above  difficulty  the  value  of  dx  in  the  present  pro- 
gram is  assumed  to  be  constant  and  the  value  of  Aty  altered  whenever  required 
during  the  calculations.  The  most  convenient  method  of  changing  Aiy  is  to 
merely  double  its  value  whenever  Eqn.  93  indicates  that  a change  is  required. 
Thus  in  the  present  program  the  stability  test  is  performed  at  each  front, 
and  if  required,  A<y  is  doubled-all  the  calculations  being  performed  within 
the  computer. 

As  a result  of  doubling  &<y  it  is  apparent  that  certain  streamline.'i,  which 
were  calculated  prior  to  the  change  of  A , are  omitted  in  subsequent  cal- 
culations. For  example,  in  the  program  the  -lines  are  numbered  consecu- 
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tively  with  the  referonco  V- line,  which  is  retained  throughout  the  calcu- 
lations, assigned  the  number  1.  Then  as  a result  of  doubling  A W the  even 
numbered  ^-lines  are  omitted  in  the  remaining  calculations.  After  the 
doubling  process  the  W -lines  are  again  numbered  consecutively.  Thus  the 

t/j -lines  numbered  N = 3,  5,  7,  before  the  doubling  process,  become 

N = 2,  3,  4 — in  subsequent  fronts.  The  general  equation  for  tracing  an 
odd  numbered  //-line  is 


/ 

n - 


n -hi 


where  n and  v!  are  the  numbered  0"— lines  immediately  before  and  after  doub- 
ling At//,  respectively. 


5.6  Detennlnatj on  of  Initial  Grid  Size 


The  initial  grid  size,  Ax  and  AW,  may  be  established  from  the  stability 
criteria  and  the  maximum  possible  number  of  grid  points  along  the  initial 
front  (x  *»  0) . Because  of  computer  storage  limitations  the  maximum  number 
of  radial  grid  points  is  limited  to  196.  In  order  to  trace  the  c^-line 
originating  at  the  outer  periphery  of  the  nozzle  it  is  advantageous  to  set 
KSLIP  = 129.  Thus,  for  given  initial  conditions  of  the  jet,  the  initial 
value  of  may  be  calculated  from  the  expression 


INlT/AL 


(<f+l)  KSLIP 


(95) 


Once  A if/  is  established  the  value  of  Ax  may  then  be  determined  from  the 
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stability  criteria.  Thus  from  Eqn.  93 


^ (96) 

A*  ~ % ,0*  u 6^L? 

where,  initially,  the  value  of  x in  Eqn.  % may  be  taken  ae  ax,  and  the 
product  /o'on"  ie  evaluated  at  the  radiue  of  the  jet  RE.  Thia  value  of 
iix  is  sufficient  to  insure  stability  at  the  first  front. 

However,  for  eubeequent  fronta  it  ia  generally  not  aufficient,  eo  that  a 
value  somewhat  smaller  than  that  given  above  ahould  be  employed.  Fran  a 
consideration  of  machine  time  the  value  of  dx  generally  ahould  be  euch 

that 

3noo 

Ax 

A ^,'la.  value  of  ax  than  that  given  above  may  be  employed  if  required; 
however,  it  ahould  be  noted  that  machine  time  may  be  increaeed  significantly, 

5.7  Tnpnt  Procedure 

In  order  to  de.cribe  the  input  procedure  it  ia  convenient  to  refer  to 
Fig.  11  wherein  ia  preaented  the  atandard  input  form  for  the  computer  pro- 
gram. The  method  of  .electing  the  value,  for  each  parameter  on  the  input 

form  is  discussed  below. 

MMTN  (Two  Options) 

The  two  option,  available  for  KMIH  are  0 or  the  number  ..signed  to  the 
(g-line  corresponding  to  the  inner  edge  of  the  mixing  region.  The 
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purpose  of  the  latter  option  is  to  permit  the  input  of  variable  flow  pro- 
perties at  the  exit  plane  of  the  nozzle  (x  = 0) , or  at  any  x-location  if 

desired* 

NMAX  (Two  Options) 

The  tv''  options  are  NMAX  = 196  or  the  number  assigned  to  the  P'-line  cor- 
responding to  the  outer  edge  of  the  mixing  region.  The  purpose  of  the 
latter  option  is  discussed  above. 

KSLIP 

KSLIP  is  the  number  assigned  to  the  ^-lir.e  originating  at  the  outer  peri- 
phery of  the  nozzle  (r  » RE).  In  most  problems  it  is  desirabls  to  set 
KSLIP  *=  129  for  purposes  of  tracing  that  jt'-line;  i.e.,  not  losing  the 
Y -line  when  ' is  doubled.  Other  values  which  may  be  employed  for 
KSLIP  are  65,  33,  17,  9,  or  5. 

NS 

The  parameter  NS  is  employed  to  permit  a curve  fit  of  the  flow  properties 
between  the  rocket  exhaust  and  secondary  streams  at  x = 0,  thereby  removing 
the  discontinuity  of  properties  at  r = RE.  The  value  of  NS  may  be  NS  - 0 
or  any  small  integer,  e.g.,  NS  « 3 is  generally  employed. 

NCHEM  (Two  Options) 

The  two  options  are  NCHEM  = 0 or  1 corresponding  to  equilibrium  or  frozen 
chemistry,  respectively. 

NCP8 

The  value  of  NCP8  is  the  number  of  temperature  values  in  the  table  of 
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TCP8  vs  Cpg  for  the  general  species. 


NXTAB 


The  value  of  NXTAB  is  the  number  of  x- values  in  the  table  of  XTAB  vs.  PTAB. 

TMDUCT  (Two  Optional 

The  two  options  are  INDUCT  = 0 or  1 corresponding  to  the  selection  of  pres- 
sure distribution  or  duct  geometry  input,  respectively. 

TSF.T  (Two  Options! 

Th.  two  options  sre  ISET  = 0 or  1 corresponding  to  the  following  conditions; 

ISET  = 0;  if  RZER01  f 0 the  program  will  stop  when  the  inner 
edge  of  the  mixing  zone  reaches  the  JKL  W-, line. 

ISET  = 1;  the  program  will  not  stop  as  above,  but  will  shift 
the  reference  C^-line  to  the  centerline  and  con- 
tiirue  the  calculations. 


The  acdsl  step  site  may  be  determined  using  the  method  outlined  in  Section 


XMAX 

The 


axial  distance  over  which  the  calculations  are  desired  is  to 


be  expressed  in  feet, 


XPRINT 


value  of  XPRINT  should  be  selected  so  that  the  ratio  XFRINT/DX  is  an 


The  value 
integer. 
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e ZEROl 

The  initial  r-value  of  the  reference  <p  -line  may  be  selected  as  0 or  any 
positive  value  depeniing  on  whether  the  reference  tf'-line  is  desired  to 
be  located  on  the  Jet  centerline  or  the  plume  slimline,  respectively  (see 
Section  5.3) • 


The  value  of  the  nozzle  radius  is  to  be  expressed  in  feet. 


EDZERO 

The  value  of  the  initial  duct  radius  is  to  be  expressed  in  feet. 


PZERO 

The  initial  pressure  in  the  secondary  stream  is  to  be  expressed  in  PSFA. 


DPDXO 

The  initial  value  of  the  pressure  gradient  may  be  input  if  known.  Other- 
wise, the  value  0 may  be  employed. 


DELTA 


The  two  options  are  0 to  1 corresponding  to  a two-dimensional  or  axisym- 
metric  flow  model,  respectively. 

TESMAX  and  TESMIN 

The  values  employed  for  TEXMAX  and  TESMIN  determine  the  width  of  the  mixing 
zone  defined  by  species  concentration  - smaller  values  yielding  wider  mix- 
ing zones.  A value  of  0.00001  is  generally  satisfactory  for  both  parameters 
(see  Section  5.4). 


REV  LTR 

U3  4288-2000  REV.  S/S4 


aa^f/s/c 


no.  D2-36251-1 


a,  FT.,  F.PSON.  and  FCON 

A complete  discussion  of  these  parameters  is  presented  in  Section  3.  Values 
generally  employed  for  the  above  parameters  are 

A “ FL  *=  1 
FCON  = 0 

The  value  of  EFSOH  must  be  determined  for  each  problem  from  empirical  ex- 
pressions  as  discussed  in  Section  3*3. 

PRT.  FLFTG.  and  FLETP 

The  turbulent  Prandtl  and  lewis  numbers  of  the  gas  may  be  specified  as  0.8 
and  1.2,  respectively.  A particle  U*.  number  of  0.6  may  be  employed  un- 
til data  become  available  for  more  accurately  establishing  its  value. 

XTAB.  RTAB.  and  PTAB 

The  data  for  the  plume  sliplin.  coordinates  and  the  pressure  distribution 
(or  duct  radius)  is  input  in  tabular  torn.  For  cases  where  the  sliplin. 
geometry  is  unknown  a value  of  RTAB  equal  to  the  nomsle  exit  radius  RE, 
should  be  employed  for  all  XTAB  values.  It  should  be  noted  that  the  input 
of  pressure  or  duct  radius  at  various  axial  distance  XTAB,  must  be  compatible 
With  the  input  of  the  control  parameter  D©UCT. 


GENU  and  Wg 

If  the  U3er  desires  to  incorporate  an  additional  chemical  species,  other 
than  any  of  those  specified  in  the  program,  the  chemical  symbol  of  the 
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component  is  input  for  GENN  and  the  corresponding  molecular  weight  for  Wg. 


TCP8  and  CP8 

If  an  additional  chemical  species  is  included  the  specific  heat  of  the  com- 
ponent as  a function  of  temperature  is  input  in  tabular  form.  It  should 
be  noted  that  this  table,  as  well  as  GENN  and  Wg,  are  input  only  if  an  addi- 
tional species  is  desired.  Otherwise,  these  parameters  may  be  omitted. 

G±,  T,  and  U (Two  Options) 

The  two  options  available  to  the  user  correspond  to  whether  or  not  it  is 
desired  to  input  variable  properties  as  the  initial  conditions.  If  the 
variable  property  option  is  selected  MEN  must  be  different  than  zero.  If 
the  constant  property  option  is  selected  the  value  of  NMIN  must  be  zero. 

The  variable  properties  are  input  as  a function  of  N.  The  constant  proper- 
ties are  input  for  the  rocket  exhaust  (C^,  T^  i^,)  and  the  secondary 

stream  TNMAX*  UNMAX^* 
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7.0  APPENDICES 

7.1  Computer  Program  Notation 

A program  notation  for  the  exponent  n,  in  Eqn.  94 

q outer  boundary  condition  for  the  concentration  of  species  i 

iNMAX 

inner  boundary  condition  for  the  concentration  of  species  i 

iNMIn 

CP8  specific  heat  of  the  general  species  (i  “ 8) 

DELTA  cf  = 0 for  two-dimensional  model 

S *=  ' 1 for  ajdsymmetric  model 

DPDXO  initial  pressure  gradient  dp/dx 

DX  axial  step  size  Ax 

EPSON  Reichardt's  constant  60,  in  Eqn.  94 

FCON  constant  C,  in  Eqn.  94 

FL  reference  length  L,  in  Eqn.  94 

FLETG  turbulent  Lewis  number  of  the  gas 

FLETP  turbulent  Lewis  number  of  the  particles 

GENN  molecular  formula  for  the  general  species  (i  =»  8) 

INDUCT  INDUCT  *=  0 for  pressure  distribution  input; 

INDUCT  = 1 for  duct  profile  input 

ISET  ISET  = 0 stops  program  when  NMIN  = 1 and  rQ  0 

ISET  — 1 program  continues  after  shifting  rQ  to  the  centerline 

KSLIP  number  assigned  to  the  ip- line  originating  at  the  lip  of 

the  nozzle  (r  = r ) 
e 

NCHEM  NCHEK  = 0 corresponds  to  equilibrium  flow; 

NCHEM  = 1 corresponds  to  frozen  flow 
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NCP8 


number  of  temperature  values  in  the  table  of  CP8  vs  T 
MiAX  Tn<nri  rmim  number  of  ip  -lines  desired  (NMAX  = 196) 

NMIN  NMXN  = 0 for  constant  property  input; 

NMIN  0 for  variable  property  input 
NS  number  of  -lines  within  which  the  flow  properties  between 

the  inner  and  outer  streams  are  curve— fitted 
NXTA£  number  of  values  in  the  table  of  X vs  PTAB 

FRT  turbulent  Prandtl  number 

PTAB  value  of  r^  (or  P)  in  XTAB  table 

PZEP.0  initial  static  pressure  in  the  outer  stream 
RDZERO  initial  radius  of  the  duct 

HE  radius  of  the  nozzle 

RTAB  radius  of  the  slipline  in  the  XTAB  table 

RZER01  initial  radius  of  the  reference  £/-line 

T static  temperature 

T temperature  values  in  the  table  of  T vs  CP8 

TESMAX  parameter  used  in  determining  NMAX 

TESMIN  parameter  used  in  determining  NMIN 

V axial  velocity 

W molecular  weight  of  the  general  species  (i  = 8) 

b 

XMAX  maximum  axial  distance  over  which  the  calculations  are 


performed 

XFRINT  axial  increments  for  print-out  of  data 
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7.2  Comparison  of  Theory  aad  Experiment 


theoretical  results  attained  frem  the  colter  program  are  capered  nerein 
vith  available  experimental  data,  to  establish!,*  *.  correlations  the 
model  employed  tor  the  eddy  viscosit,,  as  described  in  Section  4,  vas 
employed,  in  certain  cases  the  turbulent  transport  coemcients  sere  varied 
about  their  predicted  values  in  order  to  obtain  better  correlations.  In 
such  instances  the  variation  iron  the  predicted  values  is  discussed.  It 
ls  noted  that  analytical  results  were  cohered  with  data  obtained 
air  -augmentation  experiments;  hovever,  due  to  the  classified  na^re  of  the 
experimental  data  the  correlate  is  not  presented  herein  (see  Hef.  28). 

pa  the  experiments  conducted  by  ,aydev  sad  Peed  <«)  radial  static  and  total 
pressure  surveys  were  made  at  various  axial  locations  downstream  of  an  air 
Jet  exhausting  into  a quiescent  atoosphere.  Shorn  in  Figure  12  are  the 
theoretically  **  experimentally  determined  velocity  profiles  at  U.S  inches 

downs tream  of  an  air  Jet  having  a Mach  n»ber  of  1.96.  - parameter  y-  is 

, my-.  tv,e  -point  vtere  the  velocity  has  decayed 
the  radial  distance  measured  fran  the  po 

A i n+  +hP  inner  edge  of  the  nixing  zone,  cr  is  the 

to  one-half  its  value  at  the  inner  euge  u 

similarity  parameter,  *e  1=  the  effective  length  of  the  mixing  zone,  and 
y,  IS  the  velocity  at  the  inner  edge  of  the  mixing  zone.  B.  -lid  curve 
represents  the  theoretical  results  obtained  *«.  the  present  analysis  e^loy 
^ a value  of  Reiobardfs  coefficients  corresponding  to  the  value  of  - 
determined  experimentally  by  Meydev  and  Heed.  4s  cen  be  seen  the  predicted 
velocities  agree  quite  veil  with  the  experimental  values  over  the  entire 
mixing  region . Similar  element  was  obtained  between  the  experimental 
and  theoretical  velocity  profiles  at  axial  locations  nearer  the  nozzle  exit 
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In  the  experiments  of  For  stall  and  Shapiro  (19)  the  mixing  region  between 
a circular  jet  end  an  annular  coaxial  stream  was  investigated.  By  using 
10  percent  by  volume  of  helium  as  a tracer  in  the  jet,  the  profiles  of 
concentration  and  velocity  in  the  mixxng  region  were  determined.  Shown 
in  Figure  13  is  a comparison  of  the  theoretical  and  experimental  velocity 
wnd  concentration  profiles  at  an  axled,  distance  of  X/D  = 8.  As  can  be  seen 
reasonable  agreement  exists  between  the  theoretical  and  exper;  uental 
profiles.  As  noted  in  Ref.  19  the  Prandtl  and  Schmidt  numbers  were  both 
determined  to  be  approximately  0.7 > which  corresponds  to  a I^wis  number  of 
unity.  Theoretical  results  for  two  values  of  the  Lewis  number  are  shown 
for  the  concentration  profiles  in  Figure  13*  Snail  variations  of  the  lewis 
number  are  seen  to  effect  little  change  in  the  profiles. 

The  experimental  investigation  of  Zakkay,  et  al.,  was  conducted  to  determine 
the  turbulent  mixing  of  coaxial  jets  for  both  subsonic  and  supersonic  flow. 
In  establishing  a correlation,  the  data  obtained  with  argon  as  the  central 
stream  was  employed.  The  exit  Mach  number  of  the  argon  jet  was  approxi- 
mately unity.  The  Mach  number  of  the  annular  air  stream  was  1.6.  Shown 
in  Figure  l4  is  the  comparison  of  the  theoretical  and  measured  velocity 
and  concentration  profiles  at  various  axial  locations  along  the  jet  center- 
line.  In  establishing  the  correlation  the  value  employed  for  Reichardt's 
coefficient  was  £0  = 0.5  ft^/sec  from  X = 0 - 0.42  ft,  and  <Ea  - 1."  ft^/sei 
fru  X = 0.42  - 0.75.  As  can  be  seeD  good  agreement  exists  between  tae 
theoretical  and  measured  concentration  profiles  at  each  axial  location. 

The  velocity  profile  at  X = 0.75  ft  indicates  that  theory  predicts  slightly 
higher  velocity  ratios  than  those  measured  experimentally,  which  possibly 
may  be  attributed  to  an  incorrect  form  for  the  eddy  viscosity  expression. 
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However,  Judging  from  the  accuracy  of  the  predicted  concentration  profiles 
qnri  the  fair  agreement  of  the  velocity  profiles,  it  is  concluded  that  the 
present  eddy  viscosity  model  is  adequate  for  most  engineering  problems. 
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7.3  Sample  Problem 


In  order  to  illustrate  the  utility  of  the  computer  program  the  solution 
of  a typical  air-augmentation  problem  is  presented.  Consider  a rocket 
motor  with  the  following  operating  conditions: 


propellant 

solid  {16%  Al) 

nozzle  exit  radius 

0.63  ft. 

exit  velocity 

5570  ft/sec 

exit  static  pressure 

- 4 5,000.0  psfa 

exit  static  temperature 

4385°R 

density 

0.035  lbm/ft3 

Mach  number 

1.5 

'''he  secondary  stream  is  assumed  to  be  air  with  the  following  flow  condi- 
tions at  the  exit  plane  of  the  nozzle: 

static  temperature  - 1094°R 

static  pressure  - 11992  psfa 

velocity  — 856  ft/sec 

duct  diameter  - 2.52  ft. 


It  is  desired  to  perform  the  calculations  over  an  axial  distance  of  1. 575 
ft.  The  duct  is  assumed  to  be  constant  area.  The  preparation  of  the 
input  data,  presented  in  Fig.  15,  is  as  follows. 

NMIN 


Since  the  flow  properties  of  the  jet  are  assumed  constant  at  the  nozzle 
exit  plane,  the  value  of  NMIN  is  zero. 
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NMAX  is  set  equal  to  196. 


KSLIP 

The  value  of  KSLIP  is  taken  to  be 

KSLIP  = 129 


NS 

The  value  of  NS  is  arbitrarily  taken  to  be 

NS  - 3 

NCHEM 

It  is  desired  to  perform  the  calculations  for  equilibrium  flow.  Thus 

NCHEM  = 0 

NCP8 

Since  it  is  not  required  to  employ  an  additional  species,  other  than 
those  specified  in  the  program,  the  value  of  NCP8  is 

NCP8  » 0 


NXTAB 

The  number  of  x-values  in  the  XTAB  table  is 

NXTAB  = 6 


INDUCT 

It  is  desired  to  input  the  duct  geometry,  thus 

INDUCT  - 1 


1 
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beyond  the  point  vhere  the  mixing  zone 


It  is  desired  to  continue  the  run 
crosses  the  reference  streamline,  thus 

I SET  = 1 


DX 

The  value  of  DX  is  established  from  the  stability  criteria  In  the  following 
Banner.  Prom  Eqn.  95  and  96  the  initial  value  of  4 X required  for 


stability  is 


X * 


RE  ,, 

z KStiP  (.2.  LC/Pr^ 


z. 


Z ■ 12.9  (2  • A-9S 


0-  047  5 -ft 


The  value  of  jX  is  taken  somewhat  smaller  than  the  above  value  in  order 
to  maintain  a sufficiently  large  number  of  grid  points  in  the  mixing 
region  (see  Section  5*5)*  Thus 

DX  = 0.005  ft. 


X1UX 

The  value  of  XMAX  is 

XMAX  = 1.575  ft. 


XPRINT 

The  desired  Increment  between  printout  of  data  is 
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XFRINT 


0.25  ft. 


RZEP.01 

Since  the  static  pressure  of  the  secondary  stream  is  considerably  less 
than  that  of  the  rocket  exhaust  gases  it  is  desired  that  the  reference 
streamline  be  located  along  the  plume  slipline.  Thus 

RZER01  = 0.001.  ft. 

RE 

The  radius  of  the  nozzle  at  the  exit  is 

RE  - C.63  ft. 

RDZERO 

The  initial  value  of  the  duct  radius  is 

RDZERO  = 2.52  ft. 


PZERO 

The  initial  pressure  of  the  secondary  stream  is 

PZERO  = 11992.0  psfa 


DPDXO 

The  initial  pressure  gradient  is  obtained  from  the  method  of  characteristics 
program  employed  to  determine  the  slipline  coordinates  (5).  Thus 

DH)X0  = -800.0  psfa/ft, 

DELTA 

The  flow  is  axisymoetric,  thus 

DELTA  = 1.0 
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Die  4 separate  sheet  to  .Input  the  species  concentration  for 
the  first  eight  species  -with  N running  from  NMIN  to  UMAX  • 


Use  ^ separate  sheet  to  Input  the  last  two  species  concentration  and 
the  temperature  and  the  X-component  of  velocity  for  N running  from 
NMIN  to  NHAX. 
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Air  Augmented  Rocket  Mixing  Program 
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Figure  15  - Input  Form  for  Sample 

Problem 
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I 1 2 SPONSORING  MILITARY  ACTIVITY 


A theoretical  analysis  was  developed  for  predicting  the  flow  properties 
in  the  mixing  region  between  a particle-laden,  turbulent  rocket  exhaust 
and  a surrounding  air  stream.  It  was  assumed  that  the  turbulent  Doundary 
layer  equations,  modified  to  account  for  particles  were  valid  within  the 
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the  flow  were  limited  to  the  following:  H?0,  H.,  Op,  CO,  C02, _M9,  HC1, 


Al,  Alo0_ , and  one  additional  inert  species.  TKe  phenomenological  model 
employed'ior  the  turbulent  transport  coefficients  is  discussed  in  detail. 


The  solution  of  the  partial  differential  equations  was  obtained  using  the 
von  Mises  transformation,  expressing  the  equations  in  finite-difference 
form,  and  solving  the  resulting  equations  utilizing  a computer  program 
developed  for  the  SRU  1107. 


Results  from  the  computer  program  were  successfully  compared  with  experi- 
mental results  obtained  from  air-augmentation,  free  jet,  and  fuel  injection 
experiments. 
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